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One-electron 3+1 and 2+1 Dirac equations are used to calculate the motion of a relativistic elec- 
tron in a vacuum in the presence of an external magnetic field. First, calculations are carried on an 
operator level and exact analytical results are obtained for the electron trajectories which contain 
both intraband frequency components, identified as the cyclotron motion, as well as interband fre- 
quency components, identified as the trembling motion (Zitterbewegung, ZB). Next, time-dependent 
Heisenberg operators are used for the same problem to compute average values of electron position 
and velocity employing Gaussian wave packets. It is shown that the presence of a magnetic field and 
the resulting quantization of the energy spectrum has pronounced effects on the electron Zitterbewe- 
gung: it introduces intraband frequency components into the motion, influences all the frequencies 
and makes the motion stationary (not decaying in time) in case of the 2+1 Dirac equation. Finally, 
simulations of the 2+1 Dirac equation and the resulting electron ZB in the presence of a magnetic 
field are proposed and described employing trapped ions and laser excitations. Using simulation 
parameters achieved in recent experiments of Gerritsma and coworkers we show that the effects of 
the simulated magnetic field on ZB are considerable and can certainly be observed. 

PACS numbers: 31.30.J-, 03.65.Pm, 41.20.-q 



I. INTRODUCTION 



The phenomenon of Zitterbewegung (ZB) for free rel- 
ativistic electrons in a vacuum goes back to the work 
of Schrodinger, who showed in 1930 that, due to a non- 
commutativity of the velocity operators with the Dirac 
Hamiltonian, relativistic electrons experience a trembling 
motion in absence of external fields [l|. The ZB is a 
strictly quantum phenomenon as it goes beyond New- 
ton's first law of classical motion. Since the Schrodinger 
prediction the subject of ZB was treated by very many 
theoretical papers. It was recognized that ZB is due to 
an interference of electron states with positive and nega- 
tive electron energies 0, |1] . The frequency of ZB oscilla- 
tions predicted by Schrodinger is very high, correspond- 
ing to tuuz — 2mc 2 , and its amplitude is very small, 
being around the Compton wavelength h/mc — 3.86 x 
10~ 3 A. Thus, it is impossible to observe this effect in its 
original form with the currently available experimental 
means. In fact, even the principal observability of ZB in 
a vacuum was often questioned in the literature (3, 
However, in a very recent paper Gerritsma et al. [6j sim- 
ulated the 1+1 Dirac equation (DE) and the resulting 
Zitterbewegung with the use of trapped Ca ions excited 
by appropriate laser beams. The remarkable advantage 
of this method is that one can simulate the basic param- 
eters of DE, i.e. mc 2 and c, and give them desired values. 
This results in a much lower ZB frequency and a much 
larger ZB amplitude. The simulated values were in fact 



'Electronic address: Tomasz.Rusin@centertcl.pl 



experimentally observed. 

The general purpose of our work is concerned with the 
electron ZB in the presence of a magnetic field. The pres- 
ence of a constant magnetic field does not cause electron 
transitions between negative and positive electron ener- 
gies. On the other hand, it quantizes the energy spectrum 
into Landau levels which brings qualitatively new fea- 
tures into the ZB. Our work has three objectives. First, 
we calculate the Zitterbewegung of relativistic electrons 
in a vacuum in the presence of an external magnetic field 
at the operator level. We obtain exact analytical for- 
mulas for this problem. Second, we calculate average 
values describing ZB of an electron prepared in the form 
of a Gaussian wave packet. These average values can 
be directly related to possible observations. However, 
as mentioned above and confirmed by our calculations, 
the corresponding frequencies and amplitudes of ZB in 
a vacuum are not accessible experimentally at present. 
For this reason, and this is our third objective, we pro- 
pose and describe simulations of ZB in the presence of a 
magnetic field with the use of trapped ions. We do this 
keeping in mind the recent experiments reported by Ger- 
ritsma et al. We show that, employing the simulation 
parameters of Ref. @ , one should be able to observe the 
magnetic effects in ZB. 

The problem of ZB in a magnetic field was treated 
before Q, but the results were limited to the opera- 
tor level and suffered from various deficiencies which we 
mention in Appendix [F] A similar problem was treated 
in Ref. Q at the operator level in weak magnetic field 
limit. Bermudez et al. Q treated a related problem 
of mesoscopic superposition states in relativistic Landau 
levels. We come back to this work in the Discussion. 
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Our treatment aims to calculate directly the observable 
Zitterbewegung effects. Preliminary results of our work 
were published in Ref. (Toj |. 

An important aspect of ZB, which was not considered 
in the pioneering work of Schrodinger and most of the 
papers that followed it, is an existence of the 'Fermi sea' 
of electrons filling negative energy states. This feature 
can seriously affect the phenomenon of ZB, see We 
emphasize that both our calculations as well as the sim- 
ulations using trapped ions @ are based on the 'empty 
Dirac equation' for which ZB certainly exists. We come 
back to this problem in the Discussion. 

Our paper is organized in the following way. In Sec- 
tion H we use the 3+1 Dirac equation to derive the time 
dependence of operators describing motion of relativis- 
tic electrons in a vacuum in the presence of a magnetic 
field. Intraband frequency components (the cyclotron 
motion) are distinguished from interband frequency com- 
ponents (the trembling motion). In Sections III and 
IV we treat the same subject calculating averages of 
the time-dependent Heisenberg operators with the use of 
Gaussian wave packets. This formulation is more closely 
related to possible experiments. In Section V we simulate 
the 2+1 Dirac equation and the resulting electron Zitter- 
bewegung employing trapped ions and laser excitations 
in connection with the recent experimental simulation of 
electron ZB in absence of magnetic field. In Section VI 
we discuss our results. The paper is concluded by a sum- 
mary. In appendices we discuss some technical aspects 
of the calculations and the relation of our work to that 
of other authors. 



II. ZITTERBEWEGUNG: OPERATOR FORM 

We consider a relativistic electron in a magnetic field. 
Its Hamiltonian is 

H = ca x n x + ca y TTy + ca z n z + (3mc 2 , (1) 

where tt = p — qA is the generalized momentum, q is the 
electron charge, oti and /3 are Dirac matrices in the stan- 
dard notation. Taking the magnetic field B\\z we choose 
the vector potential A = (—By, 0,0). For an electron 
there is q = — e with e > 0. One can look for solutions 
in the form 

*(r) = e iKx+iKz $(y), (2) 

and we obtain an effective Hamiltonian Ti 

% = ch [(k x — eBy/K)a x + (d/idy)a y + k z a z ] + f3mc 2 . 

(3) 

Introducing the magnetic radius L = yfh/eB and £ = 
y/L — k x L we have y — £L + k x L 2 , eB/h = l/L 2 , and 
d/dy = (l/L)<9/<9£. Defining the standard raising and 
lowering operators for the harmonic oscillator 



a = (Z + d/dZ)/y/2, 
a+ = (£-3/0OA/2, 



(4) 



one has [d, d + ] = 1 and £ = (d + d + )/v / 2. The Hamilto- 
nian % reads 

mc 2 i H + E z g z 
H + E z a z -mc 2 l 



u 



(5) 



where 1 is the 2x2 identity matrix, E z = chk z , and 



H 



-Hid 



d 

d+ 



(6) 



with u> — V2c/L. The frequency to (which should not 
be confused with the cyclotron frequency lo c = eB/m) is 
often used in our considerations. 

Now we introduce an important four-component oper- 
ator 

A = di&g(A,A), (7) 

where A = diag(d, d). Its adjoint operator is 

A + = diag(i+,i+), (8) 

where A + = diag(d + ,d + ). Next we define the four- 
component position operators 

L (A ■ A). 



y 



V2 



X = — [A- A* 

iV2 



(9) 
(10) 



in analogy to the position operators y and x, see Ap- 
pendix [A] We intend to calculate the time dependence 
of A and A + and then the time dependence of y and X . 

To find the dynamics of A we calculate the first and 
second time derivatives of A using the equation of mo- 
tion: At = dA/dt = (i/h)[H,A}. Since 1 and a z com- 
mute with A and A + , we obtain 



At 



h 

i 

1 Ti \ ii. A+] 



[H, A 
H, A] 

[H,A 








(11) 

(12) 



There is (i/h)[H,A] = A t = 



1. j 

' o ' 11 consec l uence 
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and Aj = 



A t = 
At = 



A t 

A t o 

o At 
At o 



(13) 
(14) 



The second time derivatives of A and A + are calcu- 
lated following the trick proposed by Schrodinger. We 
use two versions of this trick 

Xt= {i/n)[H,A t ]= jHA t -~{H,A t }, (15) 
At = m)[H,At] - -jAfH +'-{%, At}. (16) 
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The anticommutator of At and % is 



Simple manipulations give 



Similarly 

H {n ^ } -H[ o {At,H})- (18) 

We need to know the anticommutators {H,A t } and 
{H,A+}. There is (i/h){H,A t } = lj 2 A and 
(i/h){H,A+} = -u} 2 A+, so that 



:{H,At} = 



A 
A 

2 i A+ 

o A+ 



(19) 
(20) 



Thus we finally obtain from Eqs. (|X9|) and (J20J) second 
order equations for A and A + 



A tt = (2i/h)HA t -Lu 2 A 



A^ 



-{2i/h)Atn -A+uj 2 



(21) 
(22) 



To solve the above equations we eliminate the terms 
with the first derivative using the substitutions A = 
ex-p(+iHt/h)B and A + = B + exp(— iHt/K), which gives 



Finally 



B tt = -{l/h 2 )H 2 B~iu 2 B,. 
B+ = -{l/h 2 )B + H 2 -B+uj 2 . 

B tt = -{h 2 +u 2 )B 1 
Bt = -B+((i 2 +uj 2 ), 



(23) 
(24) 

(25) 
(26) 



where ft = H/h. The solutions of the above equations 
are 



B 



-iMt 



B + = Cje 



C x + e lA1t C 2 , 

, -iMt C+ e iMt 



(27) 
(28) 



where M. = +y/ O 2 + ui 2 is the positive root of M 2 = 
Cl 2 + uj 2 . The operator M. is an important quantity in 
our considerations. Both C\ and are time-independent 
operators. Coming back to A(t) and A + (t) we have 

A{t) = e liu e~ lMt Ci + e lCtt e +zMt C 2 , (29) 
A + {t) = C+e +ljCit e- l(u +Cte~ ijkt e- i{lt . (30) 

In order to find the final forms of A{t) and A + (t) one 
has to use the initial conditions. They are 

A(0) = Ci+C 2 , 

A+(0) = C++C+, 

A t (0) = i(il-M)C 1 +i(Cl + M)C 2 , 

i+(0) = -iC+{h-M)-iC+{h + M). 



Ci = ^-Uf(0) + ^M^QAiO) + ^i(O), (31) 



C 2 = 
Similarly 



- l -M^At(Q) - ^M^nA(0) + ii(0). (32) 



C+ = --Af^M- 1 



-A+^Cim- 1 



-A+(0),(33) 



c; 



-At^M- 1 - -A+mM- 1 + -i+(0). (34) 



One can see by inspection that the initial conditions 
for *4(0) and At(0) are satisfied. It is convenient to ex- 
press At in terms of A and fl using the equation of mo- 
tion iAt = ACl — (lA. Then the first and second terms 
in Eqs. (f3"Tj) and (|33[) partially cancel out and the opera- 
tor A{t) can be expressed as a sum A(t) = Ai(t) + .4 2 (t), 
where 



Ai(t) = ~e 



iSltg— iMt 



A 2 (t) 



±_ e iQt e +iMt 



A(0) + M^Atyfl 
A(0) - M"\A(0)n 



(35) 
(36) 



Similarly, one can break A + {t) = Af(t) + A%(t), where 



Af(t) = - [A+W + CiA+WM- 1 
■#(*) = \ [A+W-CiA+^m- 1 

Using Eqs. © and (ITU1) we obtain 



e +iMt e -tfM ^ 37 j 



(38) 



j>(t) = ^(Ai(t)+A 2 (t)+At(t)+At(t)), (39) 
#(t) = -^=(i 1 (t)+i 2 (f)-i+(t)-i+(t)).(40) 

The above compact equations are our final expressions 
for the time dependence of A(t) and A + {t) operators and, 
by means of Eqs. (|3"9"|) and (|4U)) . for the time dependence 
of the position operators y(t) and X(t). These equations 
are exact and, as such, they are quite fundamental for 
relativistic electrons in a magnetic field. The results are 
given in terms of operators fi and M. To finalize this 
description, one needs to specify the physical sense of 
functions of these operators appearing in Eqs. (f35j) -([40 |) . 

As we shall see below, operators Cl and Ai have the 
same eigenfunctions, so they commute. Then the prod- 
uct of two exponential functions in Eqs. (13"5"1) - (|38I) is given 
by the exponential function with the sum of two expo- 
nents. In consequence, there appear two sets of frequen- 
cies uj + and u~~ corresponding to the sum and the differ- 
ence: uj~ ~ M. — O, and uj + ~ Ai + Cl, respectively. The 
first frequencies uj~ , being of the intraband type, lead in 
the non- relativistic limit to the cyclotron frequency w c . 
The interband frequencies ui + correspond to the Zitter- 
bewegung. The electron motion is a sum of different 
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frequency components when it is averaged over a wave 
packet. In absence of a magnetic field there are no intra- 
band frequencies and only one interband frequency of the 
order of 2mc 2 /h, see 

Each of the operators A(t) or A + (t) contains both 
intraband and interband terms. One could infer from 
Eqs. (|36|) and (|38|) that the amplitudes of interband and 
intraband terms are similar. However, when the explicit 
forms of the matrix elements of A(t) and A + (t) are calcu- 
lated, it will be seen that the ZB terms are much smaller 
than the cyclotron terms, except at very high magnetic 
fields. 

The operators (l and Xi do not commute with A 
or A + . In Eq. ([36)) the operator A acts on the expo- 
nential terms from the right-hand side, while in Eq. (|38p 
the operator A + acts from the left-hand side. The proper 
order of operators is to be retained in further calculations 
involving A(t) or A + (t). 

Let us consider the operator M 2 = VL 2 +uj 2 . Let E n /h 
and |n) be the eigenvalue and eigenvector of respec- 
tively. Then 

M 2 \n) = (Cl 2 +u J 2 )\n) = ±(E 2 + h 2 u; 2 )\n). (41) 

Thus, every state |n) is also an eigenstate of the opera- 
tor M 2 with the eigenvalue X 2 = E 2 /h 2 + uj 2 . To find 
a more convenient form of A n we must find an explicit 
form of E n . To do this we choose again the Landau gauge 
A = (—By, 0,0). Then, the eigenstate |n) is character- 
ized by five quantum numbers: n, k x ,k z , e, s, where n is 
the harmonic oscillator number, k x and k z are the wave 
vectors in x and z directions, respectively, e = ±1 labels 
the positive and negative energy branches, and s = ±1 
is the spin index. In the representation of Johnson and 
Lippman [1 21 ] the state |n) is 



obtain from Eq. ((41 



N„ 



( s 1 (eE ntkz +mc 2 ) \n - 1) 

s 2 (eE n> k z + mc 2 ) \n) 

(sip z c - s 2 fkj n ) \n — l) 

\ -(sihoj n + s 2 p z c) \n) 



(42) 



where Si = (s + l)/2 and s 2 = (s — l)/2 select the 
states s = ±1, respectively. The frequency is uj n — LOy/n, 
the energy is 



E n , kz = ^{mc 2 ) 2 + (ku n ) 2 + (hk z c) 2 , (43) 

and the norm is N nekz = (2E 2 k ^ + 2emc 2 E ritkz )~ 1 / 2 . 
In this representation the energy E n ,k z does not depend 
explicitly on s. Then the eigenvalue of operator Ct is E n = 
eE n ^ kz /h. The harmonic oscillator states are 



^ik x x+ik z z 



(44) 



27rvTC, 

where H n (£) are the Hermite polynomials and C n = 
Y / 2"n!- v /7r. Using the above forms for |n) and E n ^ kz we 



M 2 \n) = -E 2 n+1 t Jn) 



(45) 



i.e. A n = X n ,k z = ±E n+ i_ k Jh. In further calculations we 
assume \ n ,k z to be positive. The operator M 2 is diag- 
onal. As follows from Eq. (|4Tj) the explicit form of M 2 
is 



M 2 = diag[di,d 2) di,d 2 ] 



(46) 



with 



h 2 di = (mc l Y + (cpi) + h z uj z + h z u z aa + , (47) 
h 2 d 2 = (mc 2 ) 2 + (cp 2 z ) + h 2 uj 2 + h 2 uj 2 a + a. (48) 

Because M 2 = £1 2 +uj 2 , eigenstates of A4 2 do not depend 
on the energy branch index e. 

To calculate functions of operators ft and M. we use 
the fact that, for every reasonable function / of opera- 
tors & or M 2 , there is f(Cl) = J2 n f( eE n,kJfi)\n)(n\, and 
f(M 2 ) = Z n f(E 2 n+ i,kJ&)\n){ti\, see e.g 0. Thus 



< -■ = ^e ±ietB "'^ /ri |n)(n| 



±i0t 



M = (M 2 ) 1/2 = vY. 



E n +l,k z 



M- 1 = (M 



2-V-1/2 _ 



(49) 



|n)(n|, (50) 



-|n)(n|, (51) 



E n +l,k z 

e ±iMt = e ±«(^ 2 ) 1/3 = ^ e ±^ + i,^A| n )( n | ) (52) 

n 

where v = ±1. Without loss of generality we take v = 
+1. The above formulas can be used in calculating the 
matrix elements of A(t) and A + (t). 

Taking the eigenvectors |n) = |n, k x , k z , e, s) and |n') = 
\n' , k' x , k' z , e', s') with n' = n + 1, we calculate matrix el- 
ements A n ,n'(t) using A(t) given in Eqs. (|35]l and (f36"T) . 
The selection rules for ^4 n ,n'(0) are k x = k' x , and k z = k' z , 
while e, e',s,s' do not obey any selection rules. The 
matrix element of M~ 1 A(0)£l appearing in Eqs. (|35|) 
and ([35]) is 

(n|M-U(0)fi|n') = J—A(0) n . nl e ^^ = e'A(0) n , n , 

^n.k z n 

In the last equality we used E n \ kz = E n+ \, kz and X n . kz = 
E n +i,k z /h. Introducing u n , kz = E n , kz /h we obtain 
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-Al = :e«^'rV,)'(i + f ' W o) nyi (54) 

A 2 (t)n,w = l -e^^ +x ^ t (l-e')A{Q)^. (55) 

Thus the matrix element of ^4.(t) n ,n' — •Ai(t)u,a' + 
A 2 (t) n ^ n i is the sum of two terms, of which the first 
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is nonzero for e' 



-1, while the second is nonzero 



for e' = — 1 . As shown in Appendix [Bj the matrix el- 
ements obtained in Eqs. (f54)) and (|55|) are equal to the 
matrix elements of the Heisengerg operator A(t) n n i — 
(n\e mt A(0)e- mt \n'}. 

For A + (t) n > >n — Af (t)n',n + A% (*)n',n we obtain in a 
similar way 

A+(t)n>,n = ie i {- A -.».- £W -.*.) 1 (l-e / )i+(0) a> .(57) 



Formulas (f54")) -(f5"T |) describe the time evolution of the ma- 
trix elements of Ait) and -4 + (i) calculated between two 
eigenstates of Cl. The frequencies appearing in the expo- 
nents are of the form ±\ n ^ kz ±t*J n ,k z — i&)n+i,k x ±^n,k z - 
The intraband terms characterized by to c n — co n +i.k z ~ 
^n,k z correspond to the cyclotron motion, while the in- 
terband terms characterized by w% = u> n +i,k z + u n ,k z de- 
scribe ZB. Different values of e, e' in the matrix elements 
of A(t) n , n ',^2(*)n',n,^i"(*)n',n,^2 (*)n',n give contribu- 
tions either to the cyclotron or to the ZB motion. In Ap- 
pendix [B] we tabulate the above matrix elements for all 
combinations of e, e'. The exact compact results given in 
Eqs. ([5i |) -([57 ]) indicate that our choice of A(t) and A + {t) 
operators for the description of relativistic electrons in a 
magnetic field was appropriate. 

To complete the operator considerations of ZB we esti- 
mate low- field and high-field limits of A n , n i it) . Consider 
first the matrix element between two states of positive 
energies and s = —1. We take |n) = \n, k x , k z , +1, — 1) 
and In') = \n + l,k x ,k z ,+l,-l). Then 



A(t)^ nl = V^TT e^.^-^+i.fcJt/fc x 

(E n ,k z + E n+1 , kz )jE n . kz + mc 2 ) 

1\/E n , kz E n+1 , kz {E n ,k; + mc 2 )(E n+ltkz + mc 2 ) 



(58) 



This equals ^i(t) n , n ' given in Eq. (l54l) because 
y4 2 (0n,n' = for e' = +1. At low magnetic fields there is 

h 2 uj 2 HeB _ . . 

E n +i,k z ~E n ,k z = — — ~ = ftw c , (59) 

^n+i,k z + &n,k z m 

where in the denominators we approximated E n ^ kz ~ 



E n +i, k z 



and used u) = \/2c/L, see Eq. 



Set- 



ting again E n _ kz ~ E n+ \ kis ~ mc 2 in the numerator and 
denominator of Eq. (|58p we recover the well known result 
for the matrix elements of the lowering operator a in the 
non-relativistic limit 



(60) 



Consider now the above state |n) from the positive energy 
branch and the state |n') from the negative energy branch 
|n') = \n + 1, k z , k z , — 1, — !)• Then the matrix element is 



jE n ,k z — E n+ i ik JjEn tkz +mc 2 ) 

2y/E 7hkz E n+1 , kz iE n , kz + mc 2 )(E n+1 , kz - mc 2 ) 



(61) 



Assuming low magnetic fields, small k^ z values, and using 
the above approximations we obtain 



•^(^)n.n' 



hu), 



2mc 2 



c_ ~2imc 2 t/h 



(62) 



Since at low magnetic fields there is fuv c <C mc 2 , the 
amplitude of interband (Zitterbewegung) oscillations is 
much lower than that of the cyclotron motion. At low 
magnetic fields both the amplitude and the frequency of 
ZB do not depend on the quantum number n. 

Let us consider now the opposite case of very strong 
magnetic fields, when hu> 3> mc 2 and tuu hck z . Such 
a situation is difficult to realize experimentally since the 
condition haj = mc 2 corresponds to L = V2\ c , i.e. the 
magnetic length is of the order of the Compton wave- 
length. Within this limit E n ^ kz ~ E n = huj^/n, and 
the matrix elements of «4(t) n ,n' for the cyclotron and ZB 
components are 



Ait) n , n , = (y/K ± v^TT )e 4 "^ ^ 



(63) 



where the upper signs corresponds to the cyclotron and 
the lower ones to the ZB motion, respectively. 

The conclusion from the above analysis is that at low 
magnetic fields of a few tenths of Tesla the ZB amplitude 
is eight orders of magnitude smaller than the cyclotron 
amplitude. In fields of the order of 4.4 x 10 9 T the ZB 
motion and cyclotron motion are of the same orders of 
magnitude. This completes our derivation and analysis 
of the operators describing electron motion in a magnetic 
field according to the 'empty' Dirac equation. However, 
it is well known that observable quantities are given by 
average values. 



III. ZITTERBEWEGUNG: AVERAGE VALUES 

In this section we concentrate on observable quantities, 
i.e. on electron positions and velocities averaged over a 
wave packet fir). We analyze the one-electron Dirac 
equation neglecting many-body effects. Our calculations 
are first performed for a general form of f(r) and then 
specialized for the Gaussian form of the packet. 



A. Averaging procedure 



We take a packet with one or two nonzero components, 
i.e. /(r)(ai, a 2 , 0, 0) T with |ai| 2 + |a 2 | 2 = 1. According 
to the procedure adopted in the previous section, we first 
calculate the averages of A(t) and A + it) operators and 
then the position operators y(t) and X{t). We do not 
consider multi-component packets because they are diffi- 
cult to prepare and their physical sense is not clear. 

Averaging of operators A(t) and A + it) can be per- 
formed using formulas from the previous section, see 
Eqs. ([33)1 - (|3"5)) . However, a simpler and more general 



6 



method is to average the Heisenberg time-dependent 
form A(t) = e <At A(0)e _< "' with the use of two unity 
operators 1 = J2 n l n )( n l- Then the average of A(t) is 

(A(t)) = (f\A(t)\f) = (f\e i6t / h Ae- ifit / h \f) = 
^e leE '^ t/h e- ie ' E ^ t/h (n\f)(f\n')(n\A\ri'), (64) 

and similarly for (A + (t)). There is 

Yl J dk x dk' x dk z dk' z - (65) 

nn' n,n' ee' s,s' 

The selection rules for the matrix elements (n|.4|n') are: 
n' = n + 1, k' x = k x , k' z = k z , while for (n|.4 + |n') we 
have n' = n — 1, fc^, = k Xl k' z = k z . The wave packet is 
assumed to be separable /(r) = f z (z)f xy (x, y). Then we 
have 

( n l/> = Xnefc;3z( fc z)( s l a l^«-l + S2a2-Fn), (66) 

where Xnek z = (&En,fc» + mc 2 )N nekz , and 
1 



F n (k x ) 
in which 



2LC n 



gxv(Js x ,y)e-v* ft»(£)dy, (67) 



n J — oo 



and 



9xy 



9z{k z ) 



f xy (x,y)e^ x dx, (68) 



'2-k 



f z ( z y k > z dz. 



(69) 



To proceed further we must specify nonzero components 
ai,ci2 of the wave packet. First, we limit our calcu- 
lations to a one-component packet with the nonzero 
second component corresponding to the state with the 
spin s z = —1/2. Setting a\ = Q,a% = 1 we obtain from 
Eq. flU: (n|/) = s 2 Xnek^z{k z )F n {k z ). This gives 

/OO 
dk z dk' z gl{k z )g z {k' z ) x 

n,n- -°° 
y2 e i(eE nikz -e'E n , h Jt/h^ 



dk x dk' x F* n (k x )F n , (k' x ) Y S2s' 2 (n|i|n') . (70) 



The upper indices in (A(t)) 2 ' 2 indicate the second 
nonzero component of the wave packet involved. The ma- 
trix element (n|„4|n') 2,2 has ten nonzero terms. The sum- 
mation 53 s 's s 2S2(n|„4|n') gives only three nonzero terms 
being the products of (s2s' 2 ) 2 , since S1S2 = Sis' 2 = 0. 



Rearranging summations and integrations we obtain 

(A(t)) 2 ' 2 =Y,Un,n+lV^TT X 



(eE nikz -e'E rl+likz )t/h 



dk z \g z (k z )\ 2 Y^ 



[xlekzXl+Ukz + Vnek z Vn+l,ek z (<?pl + &&*)] , (71) 

where rj nekz = Xnek z N n<ikz . We define 

/oo 
F* m {k x )F n {k x )dk x . (72) 
-OO 

Since X 2 ntkz = (l/2) + emc 2 /{2E n , kz ), r lntkz = e/(2E n , K ), 
and E 2 , ^ = (mc 2 ) 2 + {cp z ) 2 + (hoj n ) 2 , we have 



(At)) 



2,2 



U n , n+ iVn+ 1 



dk z \g z (k z )\ 2 J2 e*^"-** - e ' E *+i.>°* W h x 



l + ee' 



E 



n,k z 



E n +l,k z 



E n ,k z E n +l,k z 



.(73) 



The summations over e and e' lead to combinations of 
sine and cosine functions. The calculation of (A + (t)) is 
similar to that shown above, but the selection rules for 
(n|yV~|n') are n' = n— 1, k' x = k x , k' z = k z . Performing 
the summations we finally obtain 



(At)) 



2,2 



i Y Vn+ 1 [/„,«+! 



u+ - u: 



where 
/ 



l ± 



ill + HZ 



|5 Z (^)| 2 x 



(74) 



(75) 



cos [{E n+1)kz =F E ni kjt/H] dk z , 



(76) 



= mc 2 



± 



5 Z (^)| 2 x 



E n ,k z E n+ i tkz 
sin [(£ n+ l,fc e T E nMz )t/h] dk z . (77) 

Finally, average electron positions (y(t)) 2,2 and (,£(i)) 2 ' 2 
for the 3+1 Dirac equation in a vacuum are [see Eqs. (|74[) 
and {75]), and Eqs. <[39j and @J] 

w)) 2 ' 2 = Ae^ >< 

* n 

(C/ n ,„ +1 + £/ n+1 ,„) (J+ + /") + k 0x L 2 , (78) 

* n 

+ f/n+l,n) (^ + + ^) • (79) 
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For a packet with the first nonzero component we ob- 
tain similar results. In both cases there appear the same 
frequencies but they enter to the motion with different 
amplitudes. This is illustrated in Fig. [7] of Section V for 
the 2+1 Dirac equation. The averages (y(i)) and (X{t)) 
are equal, up to a constant j/o = ko x L 2 , to the averages 
of the usual position operators (y(t)) and (x(t)), see Ap- 
pendix [AJ 

Finally we consider a two-component wave packet 
(r|/> = /(r)( ai ,a 2 ,0,0) T with K| 2 + \a 2 \ 2 = 1. Defin- 
ing fx = ax/ and f 2 = a 2 f we have 



(A + f2\A(t)\h + fr) 



= (h\Mt)\h) 
(/Mm*) - 



+ (f2\A(t)\f 2 ) 

</2li(*)l/l>, 



(80) 



and similarly for (/i + fo\A + (t)\fi + f 2 ). The first two 
terms were calculated above. The other two terms are 

(At)) 2A = \<>Zai u n.n (J+ + Jc) , (81) 

n 

(A+it)) 1 ' 2 = Ua^Un.n + J-*) , (82) 



and (Ait)) 1 - 2 = (A+{t)) 



2,1 



0. We define 



± 



Cp z hu! 



-g* z (k z )g z (k z ) x 



cos [(E n+ltkz + E nM )t/h] dk z . (83) 



The integrals Jf- describe mixing of the states with dif- 
ferent components s z . Since are odd functions of k z , 
they vanish for the wave packet with k$ z = 0. Contribu- 
tions from these integrals are relevant only for magnetic 
fields of the order of B ~ 5 x 10 9 T, where the magnetic 
length L is comparable to A c . The velocity of the packet 
in the z direction v z — hko z /m must be comparable to c. 
At low magnetic fields the mixing terms are negligible. 

All the above results were obtained the for the 3+1 
Dirac equation. A reduction to the 2+1 DE is obtained 
by setting \g{k z )\ 2 = S(k z ) in Eqs. flTJJl and ([77]) and 
performing integrations over k z . Below we quote final 
results for 0>(i)) 2 ' 2 and (X(t)) 2 > 2 for the latter case 



(3>W) 2 ' 2 
(X(t)) 2 - 2 



L 
2~3 



^ Vn + 1 (U n , n+ x + U n +i tn ) < 

n ^ 

-^j= ^2 Vn+l (U n , n+ i + U n +i, n ) | 



n+l 



mc 
E n 



mc 



E 



n+l 



sin(w^) + 



E, 
mc 



n+l 
2 



cos(<i) } +k Qx L 2 ,{8A) 



E n E. 



n+l 



sinKfi) . (85) 



r 



In the above equations we used notation E n = E n ^ z =o, 
Lo c n = (E n+1 - E n )/h and w£ = (E n+1 + E n )/h. For the 
2+1 Dirac equation the final expressions for (y{t)) 2 ' 2 and 
{X(t)) 2 - 2 are given in form of infinite sums, while for the 
3+1 DE they are given by infinite sums and integrals 
over k z . As is known from the Riemann-Lesbegues the- 
orem (see Ref. [Hj), the k z integrals over rapidly oscil- 
lating functions of time, appearing in Eqs. ()76[) and (|77j> . 
decay to zero after sufficiently long times. Therefore, the 
packet motion for the 3+1 Dirac equation has a tran- 
sient character, while that for the 2+1 DE is persistent. 
Transient and persistent ZB motions in the two cases are 
illustrated in Fig. [5] of Section V. 



hk = h(k Ox ,0,k Oz ) 

/(r) = v^m exp ("4 " 5? - i + lkor ) ■ 

The wave packet is multiplied by a four-component Dirac 
spinor (oi, a 2 , 0, 0) T . Using the definitions of g xy (k x ,y), 
F n (k x ) and U m , n , we obtain (see Refs. [l5l [lr|) 



9xy{kx ' y) = ^ e ~ idlikx ~ k0xf<r ^ (87) 

and 



B. Gaussian wave packet 



We perform specific calculations for one- or two- 



F n (k x ) = ^S^e-M^-^VM* U n (-k x c), 

y/2-KdyCn 



where D = L 2 / J L 2 + d 2 y ,c = L 3 /JL A - d$, and 



i/2 

J ^ ~ % \ 

three widths d x , d y , d z and having a nonzero momentum \J L " + V v , 



component wave packets taking the function f(r) in \P2md ( E 2 — d 2 ^ 

form of an ellipsoidal Gaussian packet characterized by A n = -j= = I ^ 2 i ^2 

(89) 
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Uu 



A* A rnj f~ -W 2 min {™,«} 



1=0 



J2 2 i n' ni V " 



x((l-(cQ) 



^ (m+ - 20/2 H m+w _ a ( -^J, (90) 



y/1 - (cQ) 




in which Q = 1/^1 +-D 2 , VF = d x DQk 0xi and F = 
d x ko x Q. For the special case of d y = L, the formula 
for Um,n is much simpler: 



TJ. — 2- 



C m C n L V 2_P 



m+n+l 



exp 



2P 2 



H 



771 



-id 2 x k Qx 
P 



, (91) 



where P = y d^ + ^L 2 . In the above expressions the 

coefficients C/ m ,„ are real numbers and they are symmet- 
ric in m, n indices. For further discussion of of U m ,n see 
Appendix [O and Ref. (l5j |. 

The coefficients t/ m .„ given in Eqs. (p?H)) and (liTTI) . apart 
from the fc z dependent parts of the integrals if and , 
describe the amplitudes of oscillation terms. In the spe- 
cial case of n — m they are the probabilities of the expan- 
sion of a packet /(r) in eigenstates of the Hamiltonian 
H = (h 2 /2m)(p — eA) 2 of an electron in a uniform mag- 
netic field. This ensures that all U n ^ n are non-negative 
and normalized to unity, so that in practice there is a fi- 
nite number of non-negligible U n ^ n coefficients. There 
is also a summation rule for \Jn + 1 C/„+i. n , see Ap- 
pendix [C] which reduces the number of non- negligible co- 
efficients U n +i tn . Finite number of non-negligible coeffi- 
cients U n ^ rn limits the number of frequencies contributing 
to the cyclotron and ZB motions. Simpler formula (j9"Tj) 
for U m ^ n shows that the coefficients U n ^ n +i are relevant if 
all the quantities d x , d y , k^ x and the magnetic length L 
are of the same order of magnitude. The remaining pa- 
rameters, i.e. d z and kg z , can be arbitrary with the 
only requirement that the total initial packet velocity 
\vq\ = h\ko\/m must be smaller than c, which is equiva- 
lent to \J k 2 x + k 2 z < A" 1 . Because of the x—y symmetry 
of our problem, it is natural to take d x ~ d y . In our cal- 
culations we keep d x d y s=s d Zl but they do not have to 
be equal. Because a constant magnetic field does not cre- 
ate electron-hole pairs, there is no restriction on B and 
the magnetic length L can be arbitrarily small. 

Before presenting numerical calculations for the mo- 
tion of a wave packet in a magnetic field we analyze 
qualitatively possible regimes of parameters for realistic 
physical situations. This problem has two characteristic 
lengths: the Compton wavelength A c = 3.86xlO _3 A and 
the magnetic length L. For a magnetic field B = 40 T 
there is L = 40. 6A. The magnetic length is equal to A c 
for B = 4.4 x 10 9 T. We then distinguish two regimes of 
parameters: i) the low-field limit, in which packet widths 
d x , d y , k^ x and the magnetic length L are of the order of 
nanometers, and ii) the relativistic regime, in which all 
quantities d x , d y , k^ x and L are of the order of A c . 



C. Low magnetic fields 

At low magnetic fields the electron moves on a circular 
orbit with the frequency uj c = eB/m and the radius r = 
mv/eB. The aim of this subsection is to retrieve the non- 
relativistic cyclotron motion from the general formulas in 
Eqs. ([75 ]) - ([75)1 . Additionally, we show that ZB exists even 
at low magnetic fields but its amplitude is much smaller 
than A c . 

At low magnetic fields we approximate E„,k z — mc 2 
and E n+ i t k z --E n ,t — Sw c - Then I~ and I~ in Eqs. ([75)1 
and ([77]) reduce to 

/oo 
\g z (k z )\ 2 dk z , (92) 
-OO 
n OO 

I-=2sm(iu c t) \g z (k z )\ 2 dk z , (93) 



and they do not depend on n. The integrals over k z give 
unity due to the normalization of the wave packet. The 
summation over n in Eqs. (|78[) - (I79I) is performed with the 
use of the formula (see Appendix [C]) 



1 



^ y/n+l U n+ i t „ = --^=k 0x L. 



(94) 



We find 



(y(t)) ~ -k 0x L 2 cos(uj c t) + k 0x L 2 , (95) 
(x(t)) ~ -k 0x L 2 sm(uj c t). (96) 

Since L 2 — h/eB and vq x = Hko x /m, we obtain ko x L 2 = 
mvo x /eB, which is equal to the radius of the cyclotron 
motion. Taking the time derivative of (y(t)) and (x(t)) 
and using definitions of L and lu c we have 



(«»(<)> 
("*(*)) 



■ sm(aj c r), 



m 
hk 0x 



cos(w c i). 



(97) 
(98) 



Thus we recover the cyclotron motion of a non-relativistic 
electron in a constant magnetic field. 

Now we turn to the ZB motion. At the low-field limit 
we again separate the integration over k z from the sum- 
mation over n. The integration over k z selects k z ~ 0, so 
the amplitude D of the ZB motion is [see Eq. (|78l) ] 



D 



2V2 
L 



1 



E n ,o 
E n -\ 



1.0 



h 2 to 2 2k 0x L 1 
x _ „ . x — = -X c (k 0x X c ) 



^ Vn + 1 (Un+l,n + U n , n+ i) 

2V2 " 2mV " V2 ~ 2" cv " us " c/ - (99) 

Thus at low magnetic fields the amplitude of the ZB mo- 
tion is a small fraction of A c , since fcoa;A c -C 1. This 
agrees with the old predictions of Lock in Ref. [l4| • An 
interesting feature of ZB motion at low magnetic fields 
is its slow decay in time, proportional to f 1 ! 2 . A simi- 
lar decay of ZB proportional to t -1 / 2 was also predicted 
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for a one-dimensional electron Zitterbewegung in carbon 
nanotubes |17]. To understand this behavior we consider 
the integral 7+ (t) in Eq. (|76[) . Retaining only the cosine 
function and taking a Gaussian wave packet we obtain 



I+(t) ~ D / cos [(£ n+ i,fc* + E n<kz )t/h] e' d ^dk z , 

J — OO 

(100) 

where 13 is a constant independent of fc z and propor- 
tional to D, as given in Eq. (|99l) . Expanding the energy 
E n ^ z in Eq. (|43l) to the lowest terms in fc z , we have 



I+(t) ~ D cos 



(2^ 

The direct integration gives 



k 2 X 



rac t 



D 



F osc (t) 



+ (ht/m) 



21 1/4' 



(101) 



(102) 



where F osc (t) is a function oscillating with the fre- 
quency uj — 2mc 2 /fi and having the amplitude of the or- 
der of unity. Therefore, the ZB oscillations decay as t^ 1 / 2 
and they persist even at times of picoseconds. This is il- 
lustrated in Fig. 2] of Section IV. 



IV. RESULTS: 3+1 DIRAC EQUATION 

We present our results for the 3+1 Dirac equation in 
a vacuum beginning with the relativistic limit for a wave 
packet with the second nonzero component. The aver- 
age packet positions y(i) and X(t), given by Eqs. (|7"8|)- 
(|79| . are calculated computing numerically the coeffi- 
cients U m ,n, see Eqs. (1UU1) and (f9"Tj) . In our calculations we 
use the first n — 400 Hermite polynomials. For each set 
of parameters L, d x ,d y , ko x we check the summation rules 
for U n ^n and yjn + 1 U n +i,n> see Appendix [Cj With the 
numerical procedures we use, these rules are fulfilled with 
the accuracy of ten or more digits. In Fig. [TJ we plot the 
electron positions calculated for the first 200i c ~ 0.25 at- 
toseconds of motion for various packet parameters. The 
time scale is in units t c = h/mc 2 = 1.29 x 10 -21 s. We 
chose magnetic field B — 4.4 x 10 9 T and an elliptic wave 
packet with fcox = 0.998A^T 1 and k$ z =0. It is seen that 
the ZB oscillations consist of several frequencies. This is 
the main effect of an external magnetic field, which quan- 
tizes both positive and negative electron energies into the 
Landau levels. At larger times the oscillations in the 3+1 
space go through decays and revivals, but finally disap- 
pear. Thus the motion of electrons shown in Fig. [TJ has 
a transient character in which several incommensurable 
frequencies appear. The calculated motion is a combi- 
nation of the intraband (cyclotron) and interband (ZB) 
components. In the relativistic regime the components 
have comparable amplitudes. The character of motion, 
number of oscillations in the indicated time interval and 
the decay times strongly depend on packet's parame- 
ters. For Fig. [TJ; we chose the packet width d y = 4.8A C . 
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FIG. 1: Calculated motion of wave packet with the sec- 
ond nonzero component during first 200 t c ~ 0.25 attosec- 
onds of motion for various wave packet parameters. The 
magnetic field corresponds to L — A c . Packet parameters: 
d x — 1.5\ c ,d z — 1.8A C , ko x = 0.998A" 1 , fco 2 = 0. Time scale 
is in t c — h/(rnc 2 ) = 1.29xl0 -21 s units, while position is in 
A c = h/(mc) = 3.86X10" 13 m units. 
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FIG. 2: Calculated motion of wave packet with the second 
nonzero component and nonzero velocity in the z direction. 
Packet parameters: d x = 2.0A C , d y — 1.8A C , d z = 1.5A C , 
k 0x = 0.673A" 1 . 
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FIG. 3: Trajectories of wave packets with the second nonzero 
component for 3+1 Dirac equation in various magnetic 
fields. Packets parameters: d x — 0. 632(5;, /B)°' 5 X C , d y — 
0.569(Bf,/B) a5 A c , d z = 0.474(B 6 /B)°' 5 A C , k 0z = 0, k 0x = 
0.999(B/B b ) - B \J l , where B b = 2 x 10 10 T. The products 
Lkox, d x kox, dykox and d z ko x are the same for all figures. In 
all cases the packet motion is transient but for lower magnetic 
fields the decay of oscillations is slow. 




0.0 2.5 5.0 7.5 t(10" 21 s) 



FIG. 4: Calculated ZB components of electron motion in a 
magnetic field in two very different time scales. Packet pa- 
rameters are: d x = 20000A C , d y = 18000A C , d z = 15000A C , 
k 0x = 0.5L" 1 = 8.72 xluV" 1 , k 0z = 0. The ZB oscillations 
decay as t _1//2 but they have very small amplitudes. Note the 
collapse-and-revival character of ZB oscillations. 



The number of oscillations is then reduced compared to 
Figs. QJi and [T}d. This confirms a previous observation 
(see Ref. [TH) that the packet parameters have to be 
carefully selected for ZB to be observable. In contrast to 
the low-field limit, in the high-field regime the amplitudes 
of ZB are of the order of A c . 

Motion of non-relativistic electrons in the z direction, 
parallel to the magnetic field, is independent of the cir- 
cular motion in x — y plane. However, the motion of 
relativistic electrons in the z direction is coupled to the 
in-plane motion. To analyze this effect we calculated 
positions of the relativistic wave packet with a nonzero 
initial velocity assuming feo = (ko x ,0,ko z ) with the con- 
straint \ko\ < X^ 1 . In Fig. [5] we show the calculated 
packet motion with fixed ko x = O.GliX^ 1 and various 
values of k 0z for B = 4.4 x 10 9 T. As seen in Figs.[Da,,[2b 
and [2k, the existence of nonzero ko z component reduces 
the number of oscillations in the cyclotron and ZB mo- 
tions. Increasing k$ z leads to a faster decay of the mo- 
tion. The maximum initial amplitudes of oscillations do 
not depend on ko z , but the amplitudes at larger times 
decrease with increasing k§ z . 

To visualize the gradual transition from the non- 
relativistic to the relativistic regime we plot in Fig. [3] 
the packet trajectories for four values of magnetic field. 
In all cases the packet parameters are chosen in a system- 
atic way keeping constant values of the products: Lko x = 
0.47, d x k 0x = 0.632, d y k 0x = 0.569 and d z k 0x = 0.474. 
For B = 2 x 10 7 T the trajectories of electron motion are 
still circular, as at low magnetic fields. When the field 
is increased to B = 2 x 10 s T the trajectories are de- 
formed into slowly decaying spirals. At very high fields: 
B = 2 x 10 9 T and B = 2 x 10 10 T, the trajectories are 
described by fast decaying spirals. The amplitude of mo- 
tion decreases with increasing field, which is caused by 
the decrease of magnetic length L. 

Finally, in Fig. @] we plot the ZB part of motion at the 
low magnetic field B = 20 T in two scales of time. The 
amplitude of ZB motion is D = 6.5 x 10~ 8 A, which agrees 
well with its estimation given in Eq. In Fig. 0^ we 

observe a slow decay of oscillations with its envelope de- 
caying as t" 1 / 2 . In Fig. |3b we show rapid ZB oscillations 
with the frequency uj z = 2mc 2 /h = 7.76x 10 20 s _1 . The 
ZB oscillations exist even at times of the order of picosec- 
onds. 

Generally, the ZB effects are observable in magnetic 
fields of the order of 4.4 x 10 9 T for wave packets mov- 
ing with an initial velocity close to c. These packets 
should have width of the order of A c . It is not possible 
to fulfill all these requirements using currently available 
experimental techniques. In addition, the predicted am- 
plitudes of the ZB motion are of the order of A c , which 
makes their experimental detection extremely difficult. 
However, there exists now a very powerful experimental 
possibility to simulate the Dirac equation and its conse- 
quences. We explore this possibility in the section below. 
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FIG. 5: Calculated motion of two-component wave packet 
simulated by trapped 40 Ca + ions for three values of effective 
rest energies HQ. Trap parameters: r\ — 0.06, Ct = 2n x 
68 kHz, A ~ 96A; packet parameters: ko x = A" 1 , d y — 
Ay/2, d x — 0.9d y . Simulations correspond toK= hcj c /2mc 2 = 
16.65 (a), 0.26 (b), 0.116 (c), respectively. Positions are given 
in L — y/2A units. Oscillations do not decay in time. 



SIMULATIONS BY TRAPPED IONS 



FIG. 6: Trajectories of electron wave packet in a constant 
magnetic field for various simulated rest energies Ml, as cal- 
culated for 2+1 Dirac equation. Trap and packet parameters 
are the same as in Fig. [5] Positions are given in magnetic ra- 
dius L. In the non-relativistic limit (a) the ZB is practically 
absent. As the rest energy decreases, the motion becomes 
more relativistic and the ZB (interband) frequency compo- 
nents become stronger. The ratio k defined in Eq. (|110l) is: 
(a) 0.0018, (b) 0.029, (c) 0.116, (d) 1.05. 



The main experimental problem in investigating the 
ZB phenomenon in an external magnetic field is the fact 
that, for free relativistic electrons in a vacuum, the ba- 
sic ZB (interband) frequency corresponds to the energy 
hujz — 1 MeV, whereas the cyclotron energy for a mag- 
netic field of 100 T is huj c ~ 0.01 eV. Thus the magnetic 
effects in ZB are very small. However, it is now possible 
to simulate the Dirac equation changing at the same time 
its basic parameters. This gives a possibility to strongly 
modify the critical ratio huj c /2mc 2 making it more ad- 
vantageous. In the following we propose how to simulate 
the 3+1 and 2+1 Dirac equations in the presence of a 
magnetic field using trapped ions and laser excitations. 

First, we transform the Dirac equation to the off- 
diagonal form 



ii = ajpi + 8mc 2 , 



(103) 



using the unitary operator P = 5(5 + /3)/v2, where 5 = 
(3 [l8[ . After the transformation the Hamiltonian 
H' 
ff't 



is ft = 



where 



H = 



cp z — irac cp x 



cp x 



hu>a 



+ 



cp z 



tUUCly 

- imc 2 



(104) 



and a y and a+ are given in Eq. (U]). 

Next we use the procedures developed earlier and con- 
sider a four-level system of Ca or Mg trapped ions [l9r - 
|2~H . Simulations of cp x and cp z terms in the above Hamil- 
tonian are carried out the same way as for free Dirac 
particles using pairs of the Jaynes-Cumminngs (JC) in- 
teractions 



a a e 



(105) 



and the anti- Jaynes-Cumminngs (AJC) interactions 



fr<t>b 
n AJC 



hrjhia+a+e 1 ^ + a~ ae" 1 ^). (106) 



A simulation of mc 2 is done by the so called carrier in- 
teraction 



H c = m{a + e^" + a 



(107) 



Here ft and fl are coupling strengths and r\ is the Lamb- 
Dicke parameter [l9j . The operators a and a + are low- 
ering and raising operators of the one-dimensional har- 
monic oscillator, respectively. These operators can be 
associated with the three normal trap frequencies and, 
therefore, with the motion along the three trap axes. Set- 
ting pairs of lasers beams in the x, y and z directions it 
is possible to simulate the lowering and raising operators 



12 




0.0 0.5 



1.0 



1.5 2.0 2^5 1 (ms) 



FIG. 7: Simulated ZB motion of one-component packets in 
the regime k> 1. Simulated gap frequency is Q, = 2-zr x 1000 
Hz, other trap and packet parameters as in Fig. [S] Upper 
part — packet with the first nonzero component; lower part 
— packet with the second nonzero component. Note largely 
different magnitudes of the x oscillations in the two cases. 




60 t (ms) 



FIG. 8: Collapse and revivals of packet motion for simulations 
using 3+1 DE (a) and 2+1 DE (b). Packet parameters: d x — 
d y — d z — L, ko x — A~ . Trap parameters as in Fig. [5] 
simulated gap frequency is f2 = 2n x 12000 Hz. Note transient 
character of motion for the 3+1 DE and persistent oscillations 
for the 2+1 DE. In both cases the collapse and revivals appear. 



along these directions, respectively. As an example of this 
procedure, one selects a pair of JC and AJC interactions 
in the x direction, adjusting their phases 4> r — —ir/2 and 
(fib = +7r/2. This way one can simulate the 2x2 Hamil- 



tonian HZ* 



H Tc + h jc t0 § et 



H% = ih v n<j x (a+ - a+) = 2r) q na x A xPx , (108) 

where p x = ih(a^ — a x )/A x . Using this technique, the p x 
and p z dependent parts of the Dirac Hamiltonian (|104[) 
can be simulated by appropriate combinations of JC and 
AJC interactions. On the other hand, a simulation of a y 
and CLy operators (which include the magnetic field) can 
be done by single JC or AJC interactions. Using the 
notation of Refs. [l^-[lll one simulates the complete 3+1 
Hamiltonian T-L by the following set of excitations 



O-x(OC) 



TT<Pr=K 

JC(ad) 



TT<Pb=K 

^AJCibc) 



H p * 



al(ac) ~ H al(bd) + H v y (ac) + H a y {bd)i ( 109 ) 



where B a ) = 2i] q QajA q p q , p q = ih(a+ - a q )/ A q , j,q = 
x, z. The subscripts in parentheses of Eq. (|109[) symbolize 
states involved in the transition in question. The spread 
of the ground ion wave function is A q = y/h/2Mv q and 

the Lamb-Dicke parameter is r] q = kyJh/2Mv q , where 
M is ion's mass, v q is trap's frequency in the q direction 
and k is the wave vector of the driving field in a trap. The 
JC interaction gives a y in H 12 and ai in H 2 \ elements of 
the Hamiltonian H in Eq. (1104[) . respectively, while AJC 

gives a y in H 2 i and a y in H^ 2 elements, respectively. A 
simulation of the 3+1 DE by Eq. (|109|) can be realized 
with 12 pairs of laser excitations: two pairs for each of the 
four interactions simulating p x and p z terms and one pair 
for each of the four remaining terms. If one omits the p z 
interaction, which corresponds to the 2+1 DE, one needs 
8 pairs of laser excitations: two pairs for the p x terms and 
one pair for the each of four remaining terms. Simulated 
magnetic field intensity can be found from the following 
correspondence (see Ref. [19(): a y — a+ = \/2L(d/dy) = 

2A(d/dy), which gives L/V2 A, where A x = A y = 
A z = A. Since the other simulations are: c ^ 2-qACl and 
mc 2 Ml, we have for the critical ratio 



heB I i]n 

m(2mc 2 ) \ 



(110) 



Therefore, by adjusting the frequencies f2 and £1 one 
simulates different values of k — hoj c /2mc 2 . This il- 
lustrates the fundamental advantage of simulations by 
trapped ions. 

In Fig. O we show the calculated Zitterbewegung 
for three values of n: 16.65, 1.05, 0.116, us- 
ing a two-component electron wave packet (r|/) = 
f(r)(V2/2,V2/2,0,0). The electron motion is a combi- 
nation of (;P) M (t), (y) 2 ' 2 {t), (y) hl (t) and (X) 2 ' 2 (t) com- 
ponents. There are no mixing terms of the form (3^) 1:2 (i) 
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etc., since they vanish for the 2+1 Dirac equation due to 
their proportionality to p z , see Eq. (1531 . The essential 
feature of the simulated characteristics is their low fre- 
quency and large amplitude of ZB. Further, it is seen 
that, as K gets larger (i.e. the field intensity increases 
or the effective gap decreases), the frequency spectrum 
of ZB becomes richer. This means that more interband 
and intraband frequencies contribute to the spectrum. 
Both types of frequencies correspond to the selection 
rules n' = n ± 1. Thus, for example, one deals with ZB 
(interband) energies between the Landau levels n = 
and nf = 1, and n — 1 and n' — 0, as the strongest 
contributions. For simulated high magnetic fields corre- 
sponding to k > 1, the interband and intraband com- 
ponents are comparable and one can legitimately talk 
about ZB. We believe that the ZB oscillations of the 
type shown in Fig. [SJi, resulting from the 2+1 DE for 
k = hio c /2mc 2 > 1, are the best candidate for an obser- 
vation of the simulated trembling motion in the presence 
of a magnetic field. The calculated spectra use the trap 
and wave packet parameters already realized experimen- 
tally, see |6| . We emphasize the tremendous differences of 
the position scales between the results for free electrons 
in a vacuum, shown in Fig. [TJ and the simulated ones 
shown in Fig. [5] The anisotropy of ZB with respect to 
(x(t)) and (y(t)) components, seen in Figs. [T] and [5l is due 
to the initial conditions, namely ko x ^ and ko y = 0. A 
similar anisotropy was predicted in the zero-gap situation 
in graphene [15j . 

In Fig. [6] we show the results of our calculations for 
different M2, simulating effective values of 2mc 2 , at a 
constant simulated magnetic field. Packet parameters 
are the same as in Fig. [5] The results are shown for ini- 
tial time intervals of the motion. In the non-relativistic 
limit illustrated in Fig.[S^, the motion is completely dom- 
inated by the intraband frequencies and it represents a 
cyclotron orbit. As the gap decreases, the motion is more 
relativistic and the circular trajectories turn into spirals. 
Simultaneously, the interband Zitterbewegung frequen- 
cies come into play. In highly relativistic regimes (low 
values of Ml) the trajectories look chaotically. However, 
the motion is not chaotic, it consists of a finite number 
of well defined but incommensurable frequencies. The 
illustrated motion of the wave packet for the 2+1 Dirac 
equation is persistent, its amplitude experiences infinite 
series of collapse and revival cycles. In the relativistic 
regime the motion is somewhat anisotropic with respect 
to the x and y directions which is related to the initial 
conditions fco = (&0x>0). This phenomenon has an anal- 
ogy in the field-free case for the relativistic regime, where 
the ZB oscillations occur in the direction perpendicular 
to the initial packet velocity pj], [H| . 

In Fig. [7] we analyze ZB of the one-component packets 
having a non- vanishing first or second component. Inter- 
estingly, they look distinctly different, and the x parts 
of the motion have different limits for mc 2 — > (i.e. for 
very small energy gaps M7). The y components of mo- 
tion are comparable in both cases, but the x components 



differ substantially. 

In all the figures presented above we showed the packet 
motion in short time spans. In Fig.[8jwe analyze the long- 
time packet evolution according to the simulated 3+1 
and 2+1 Dirac equations. In both cases the collapse and 
revival cycles occur. However, the motion according to 
the 3+1 Dirac equation is decaying in time, while the 
oscillations in the 2+1 case are persistent in time. 



VI. DISCUSSION 

We briefly summarize the important new effects 
brought to ZB by an external magnetic field: (1) The 
quantization of the spectrum for positive and negative 
electron energies results in numerous interband frequen- 
cies contributing to ZB, (2) The presence of B intro- 
duces an important new parameter into the phenomenon 
of ZB affecting all the frequencies, (3) The presence of 
intraband frequencies raises the question of what should 
be and what should not be called ZB. In our opinion, 
the interband frequencies are the signature of ZB while 
the intraband frequencies (the cyclotron resonance in our 
case) are not, (4) The presence of B 'stabilizes' ZB in the 
2+1 case making it a stationary phenomenon, not decay- 
ing in time. The last feature is related to the fact that 
the magnetic field is represented by a quadratic potential 
and, as is well known, the wave packet in a parabolic po- 
tential is not spreading in time. However, a slow decay 
of ZB in time might occur if the trembling electron emits 
radiation. This does not occur if the electron is in its 
eigenstate but it will happen if the electron is prepared 
in the form of a wave packet, because the latter con- 
tains numerous eigenstates of the electron in a magnetic 
field, see Eq. (|4"2"j). The emitted radiation can have multi- 
pole character depending on the electron energy [23T4251 ] . 
it may also be due to spontaneous radiative transitions 
between various Landau levels in the strictly quantum 
limit. Finally, in the classical limit of very high electron 
energies one may deal with the synchrotron radiation, ra- 
diative damping, etc., but this limit is beyond the scope 
of our paper. Also, a broadening of Landau levels due to 
external perturbations results in a transient character of 
ZB, c.f. (a!]. 

The time-dependent electron motion, as obtained in 
the operator form [see Eqs. (|31?)) - (l4*(Jl) ]. is described by 
four operators. We show in Appendix lEl that these oper- 
ators have different limits for low magnetic fields. How- 
ever, all of them contain both interband and intraband 
frequencies. Thus, in both operator and average formula- 
tions the cyclotron and trembling motion components are 
mixed. The method of direct averaging of operators in 
the Heisenberg form, used in Section III, is simpler than 
that of averaging the explicit forms of A and A + , as de- 
rived in Section II, since it does not require the detailed 
knowledge of these operators. The main disadvantage 
of the direct averaging is that it obscures the detailed 
structure of electron motion shown in Eqs. (|3T) )) - (|40I) . 
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In our considerations we used one-component and two- 
component wave packets and showed that the character 
of ZB oscillations in the two cases is similar, but not iden- 
tical. Calculation for three- and four-component packets, 
although possible, are much more complicated not intro- 
ducing anything new at the physical level. 

High magnetic fields for relativistic electrons in a vac- 
uum are often characterized by the so-called Schwinger 
critical field B cr for which HeB/m — mc 2 or, equiva- 
lently, L — [h/eB) 1 / 2 = h/rnc = A c . This corresponds 
to the gigantic field B cr = 4.4 x 10 9 T existing only in 
the vicinity of neutron stars. However, in simulating the 
analogous situations in semiconductors [28[ or by trapped 
ions [6[ , the corresponding critical fields are not high and 
they depend on parameters of the system in question. 
We emphasize that our results are not limited by any 
particular value of B and they describe both weak and 
high field limits. 

As mentioned in the Introduction, the initial Dirac 
equation JT]) and our resulting calculations, as well as the 
simulations based on trapped ions, represent the 'empty' 
Dirac Hamiltonian which does not take into account the 
'Fermi sea' of electrons in a vacuum having negative ener- 
gies. This one-electron model follows the original consid- 
erations of Schrodinger's. The phenomenon of electron 
ZB in a vacuum is commonly interpreted as resulting 
from an interference of electron states corresponding to 
positive and negative electron energies. The character- 
istic interband frequency of ZB is a direct consequence 
of this feature. The initial electron wave packet must 
contain these positive and negative energy components. 
It may be difficult to prepare such a packet if all nega- 
tive energies are occupied. What is more, the fully occu- 
pied negative energies mayprevent the interference (and 
hence ZB) to occur, see [a, It has been a matter 

of controversy what happens when an electron-positron 
hole pair is created by a gamma quantum [ll|. On the 
other hand, a system with negative electron energies can 
be relatively easily created in semiconductors, see [28j . 
It should be mentioned that an external magnetic field 
does not create by itself the electron-positron pairs. We 
emphasize again that our present calculations and the 
experimental simulation of Ref. [6| are realized for the 
one-electron Dirac equation for which ZB certainly ex- 
ists. 

Bermudez et al. Q treated the problem of time de- 
pendent relativistic Landau states by mapping the rel- 
ativistic model of electrons in a magnetic field onto a 
combination of the Jaynes-Cummings and Anti-Jaynes- 
Cummings interactions known from quantum optics. For 
simplicity the p z — restriction was assumed. Three 
regimes of high (macroscopic), small (microscopic) and 
intermediate (mesoscopic) Landau quantum numbers n 
were considered. In all the cases one interband frequency 
contributed to the Zitterbewegung because the authors 
did not use a gaussian wave packet to calculate average 
values. 

Our exact calculation of Zitterbewegung of relativistic 



electrons in a vacuum in the presence of a magnetic field 
and its simulation by trapped ions are in close relation 
with the proof-of-principle experiment of Gerritsma et 
al. @, who simulated the 1+1 Dirac equation and the 
resulting electron ZB in absence of magnetic field. Our 
results show that, paradoxically, the simulation of the 
DE with a magnetic field is simpler than that without 
the field. However, there is a price to pay: one needs 
at least the 2+1 DE to describe the magnetic motion 
since B parallel z couples the electron motion in x and y 
directions. 



VII. SUMMARY 

In summary, we treated the problem of electron Zit- 
terbewegung in the presence of a magnetic field in three 
ways. First, we carried calculations at the operator level 
deriving from the one-electron Dirac equation the exact 
and analytical time-dependent equations of motion for 
appropriate operators and finally for the electron trajec- 
tory. It turned out that, in the presence of a magnetic 
field, the electron motion contains both intraband and 
interband frequency components, which we identified as 
the cyclotron motion and the trembling motion (ZB), 
respectively. Next, we described the same problem us- 
ing averages of the Hcisenbcrg time-dependent operators 
over Gaussian wave packets in order to obtain physical 
quantities directly comparable to possible experimental 
verifications. We found that, in addition to the usual 
problems with the very high frequency and very small 
amplitude of electron Zitterbewegung in a vacuum, the 
effects of a magnetic field achievable in terrestrial con- 
ditions on ZB are very small. In view of this, we simu- 
lated the Dirac equation with the use of trapped atomic 
ions and laser excitations in order to achieve more favor- 
able ratios of (heB /m) /(2mc 2 ) than those achievable in 
a vacuum, in the spirit of recently realized experimental 
simulations of the 1+1 Dirac equation and the result- 
ing electron Zitterbewegung. Various characteristics of 
the relativistic electron motion were investigated and we 
found that the influence of a simulated magnetic field 
on ZB is considerable and certainly observable. It was 
shown that the 3+1 Dirac equation describes decaying 
ZB oscillations while the 2+1 Dirac equation describes 
stationary ZB oscillations. We hope that our theoreti- 
cal predictions will prompt experimental simulations of 
electron Zitterbewegung in the presence of a magnetic 
field. 



Appendix A 

In this Appendix we briefly summarize the similarities 
and differences between operators V and X ', as defined in 
Eqs. fl9|- ([T0]) . and the position operators y and x. The 
operators y — (L/%/2)(a + a + )diag(l, 1, 1, 1) and X = 
(L/iy/2)(a — a + )diag(l, 1, 1, 1) are 4x4 non-commuting 
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Operator 


(+1.+1) 


(+1,-1) 


(-++1) 


(-1,-1) 


[e int Ae- int ] n , n , 

Al(t) a , n > 
A (*)»,»' 


Mw n -u> n ,)t A 
o "mi ,11 
i(w n -w n /)t J 







e »{-%-v) f l , 

C •'Ml, 11 

t(-w n -w n ,)f J 




pi(-W n +W n /)t j 

c •^ l n,n' 



i(-B„+U„,)t J 

c «' l n.n 


[e mt A + e- m \,^ 

■^4-1 (^)ii'.n 


ir ,n 






<.*(-"»' -«»>U + , 

n',n 


n' ,n 




i(-w /+w«)< i + 



i(-u n /+u„)tj+ 

n' .n 



TABLE I: Three upper rows: matrix elements of the Heisenberg operator A(t) — e' nt A(0)e~' at and matrix elements of the 
explicit form of A(t) — Ai(t) + A2{i), as given in Eqs. (|35|) and (|36|l . calculated for four combinations of (e,e')- Three lower 
rows: the same for the operator A + (t) = e iClt A(0)e~ int and the explicit form A + {t) = At(t) + At(t). 



matrices: = 1, while the position operators y, x 

obviously commute. However, the matrix elements of y 
and X between states |n) and |n'), given in Eq. (1421) . 
are equal (up to a constant ya — k x L 2 ), to the matrix 
elements of y, x between the same states. 

As an example of this property we calculate the ma- 
trix elements of y, X, y, x at t = between two states 
|n) = \n,k x ,k z ,e, -1) and |n') = \n', k' x , k' z , e', -1) given 
in Eq. (g2]). We have 

(n|3>|n'> = -^={(n\a + a + \n')(xnXn'+N n N n ,c 2 p 2 z ) 

+ (n- l|o + a + |n' - 1) N n N n , h 2 u n u n/ } ,(Al) 

where \n) is defined in Eq. (1441) and we omitted indices 
k x and k z . For the matrix element (n|y|n') we obtain the 
same expression as in Eq. (|A1|) but with (L/y2)(a + a + ) 
replaced by y. Because (L/y/2)(a + a + ) = y — k x L 2 we 
obtain from Eq. (IA1|) 

(n|j>|n') = (n\y\n') -(n\k x L 2 \n')( XnX n'+N n N n ,c 2 P 2 z ) 
- (n - l\k x L 2 \n' - l)N n N n >h 2 uj n uj n i 
= (n\y\n')-k x L 2 , (A2) 

In order to calculate the matrix elements of x we 
observe that the Hamilton equations give: x = ca Xl 
y = ca yi p x — and p v — ca x eB. From the above 
relations one obtains p y = eBx = (h/L 2 )x, which gives 
after the integration over time 

x(t) = (L 2 /h)p y (t) + D. (A3) 

The constant of integration D can be set equal to zero 
by an appropriate choice of x(0). Since p y = (h/i)d/dy 



with d/dy = (l/L)d/d£ and d/d£, = (6 - a+)/V2 [see 
Eq. ©], there is p y {t) = (h/iLV2){A{t) - Jtf(t)), sec 
Eqs. ©-©. Thus we have 

(n\x(t)\n>) = (n|^(i(t) - A (t))|n') = <n|*(t)|n'>. 

(A4) 

Since A and are four-component lowering and raising 
operators, the selection rules for x and for X are n' = 
n ± 1 , k x — k' x and k z = k' z . There is no selection rules 
for e,e' and for s,s'. Equations (|A2[) and (IA4[) are the 
required relations between the matrix elements of y, X 
and y, x operators, respectively. 

For the states |n) and |n') with s = +1 there is also 
(ji|DV) = (n|y|n') - k x l? and (n|^|n') = (n|i|n'). For 
the states |n) and |n') with different spin indexes s and 
s' the constant term ya — k x L 2 does not appear. 

Finally we calculate the average values of y, x, y and X 
operators using a Gaussian wave packet |/) from Eq. (|86l) . 
At t = there is (f\y\f) = and (f\x\f) = 0. Next, 

(f\X\f) = L{f\^\f) = L^(f\^\f) = 0, (A5) 

and 

(f\y\f) = (f\v\f) - (f\k x L 2 \f) = -k 0x L 2 . (A6) 

All figures above refer to the averages (y(i)) and (X(t)) 
i.e., equivalently, to (y(t)} — ya, (x(t)}, respectively. 



Appendix B 



We want to prove equivalence of the general Heisenberg 
form of operators A{t) = e int A(0)e~ lQt and their explicit 
time-dependent form given in Eqs. (1331) and (|3"fJ)) . We do 
this by showing that the matrix elements of A{t) ob- 
tained by the Heisenberg formula and by using Eqs. (|35|) 



and (|36|) are the same. To calculate the matrix ele- 
ments we take two eigenstates of the operator Cl: |n) = 
| n., k x , k z , e, s) and |n') = |»i', k' x , k' z , e', s') with n! = n + l. 
We use Eq. (15^1) for the matrix element of Ai(t) and 
Eq. ([55]) for the matrix element of ^(f). On the other 
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hand, we calculate the matrix elements of e int A(0)e~ int . 
We compare the matrix elements calculated by the two 
methods for all combinations of the band indexes e, e'. 
Writing w n = E, hk JH, u n , = E n >^ k Jh, and A„ :fci = w n > 
we obtain results summarized in Table 1. It is seen that 
the matrix elements of A(t) = e int A(0)e~ int are equal 
to the matrix elements of A{t) = A\ it) + Ai (t) . Since 
the states |n) form a complete set, the equality holds for 
every matrix element of A{t). This way we proved the 
equivalence of the two forms of A(t). It is to be noted 
that selecting v — — 1 instead of v = +1 in the definition 
of the square root of operator M 2 , see Eqs. P§j) . leads 
to the same results. 

Appendix C 

Here we consider some properties of the coefficients 
U m ,n, as defined in Eq. ([72]) . First, we prove the sum 
rule J2 n U n ,n = 1- Let \n,k x ) be an eigenstate of the 
Hamiltonian H = (h 2 /2m)(p~ eA) 2 . In the standard no- 
tation there is (r\n,k x ) = e lk - x R n ^)e~^ 2 I 2 /VLC n . For 
any normalized state |/) we have 

°° />oo 

l = {f\f) = Y, dk x (f\n,k x )(n,k x \f). (CI) 

n=0 J -co 

Since F n (k x ) = (n, k x \f), see Eq. (|67|) . we obtain 

oo /-oo oo 
1 = J2 K(kx)F n (k x ) dk x = Un,n- (C2) 

n J — OO n 



n=0 ' 



n=0 



This proves the normalization of U n ^ n . Since the integral 
in Eq. (|C2j) can be expressed as /_ \F n (k x )\ 2 dk x , it 
is seen that U n , n are non-negative. The above sum rule 
was used to: i) verify the accuracy of numerical computa- 
tions of U m _m h) estimate the truncation of infinite scries 
appearing in the calculation of y(t) and X(t). 

Now we calculate another sum rule. Consider an av- 
erage value J of the operator a + over a two-dimensional 
wave packet J = (f xy \a + \f xy ) . Inserting the unity oper- 
ator 1 = J2 n I dk x \n, k x ){n, k x \ we have 

J = (fxy\d^\f X y) = ^ K / dk x (f X ya~*~\n,k x )(n,k x \f X y}. 

n=0 J-co 

(C3) 



Using the definitions of F n (k x ) and U mtn [see Eqs. (|67|) 
and (|72|) ] we obtain 

00 />oo 

y \n + 1, k x ){n, k x \f xy )yjn + 1 dk x 

n=0 J-oo 
°° roc 

= J2 / V^+T F* +1 (k x )F n (k x ) dk x 

n=0 J-co 

OO 

= ^\^TT(/„ +1 ,„. (C4) 



To calculate J independently we take the wave packet 
=^7 exp ^ — 

x y 



fxy {&> y) 



v^a exp (~^ 2 ~h +lkoxX 



(C5) 



and calculate J inserting the unity operator 1 
/ dk x \k x ){k x \. This gives 



J = 



{fxy\kx)d~^ (k x \f X y) dk x dy 



= J 9*x V ( k x,y)^ U - ^) 9xy{k x ,y) dk x dy. (C6) 

Since £ = y/L — k x L, and <9/<9£ = Ld/dy, the integrations 
over d y and fc^ are elementary and we find 



n+l.n 



ri=0 



V2 ' 



(C7) 



The above sum rule was used for an additional verifica- 
tion of [/„+i. n terms and for the analytical calculation 
of motion of a non-relativistic electron, see Eqs. ([55)1 
and (I96i 



Appendix D 

Here we calculate the average electron velocity, limit- 
ing our discussion to a packet with the second nonzero 
component. The x and y components of the velocity 
are the time derivatives of (X(t)) 2 ' 2 and (y(t)} 2 ' 2 . Since 
(X(t)} 2 > 2 and (y(t)) 2 ' 2 are combinations of (A{t)) 2 > 2 and 
{A + {t)) 2 > 2 [see Eqs. ^ and gD])] we calculate the time 
derivatives of (A(t)) and (A + (t)), as given in Eqs. ([71)1 
and ([751) . respectively. The average velocities are 



(«»(*)> 



2.2 



07+ 9/, 



<9i 



0t 



(Dl) 



L 

dl+ dl. 



22 Vn + l (U rhn+ i + U n+ i tn ) 



Of 



Of 



(D2) 



where 



die r f°° ^ „ M? 
L-^f- = ±V2c / \g z (k z )\ 2 x 

Ci J-00 Ai+l,*:* 



L 



dt 



TV2c J 



sin [{E n+1M T ^n.fc* )t/ft] dk z , (D3) 

00 2 1 

mc niu 



9z{k z )\ z x 



n=0 



cos [(E n+ i, kx T E n:k Jt/h] dk z . (D4) 

In the above equations we used E 2 +1 kz — E 2 kz — 
H 2 uj 2 . It is seen from Eqs. (|D3|) and (|Di| that the in- 
tegrals (Ldl~ /dt), describing the cyclotron motion, and 
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the integrals (Ldl+ /dt), corresponding to the ZB mo- 
tion, have the same factor (hio/E nj k z )\g z (k z )\ 2 . Inte- 
grals (Ldl~ I dt) and (Ldlf / dt) have the same property. 
Therefore the amplitudes of the cyclotron velocity and 
the ZB velocity are of the same order of magnitude. On 
the other hand, the amplitudes of positions differ by sev- 
eral orders of magnitude. 

Alternatively, we calculate the average velocities for 
the canonical velocity operators. The velocity operator 
is obtained from the equation of motion v = (i/h)[H,f], 
which gives v x = ca x and v y — ca y . Now we show 
that the average velocities obtained in Eqs. (|Dl[) - (|D2j) 
are equal to the averages of v y (t) and v x (t) . We limit our 
calculations to a wave packet with the second non-zero 
component. 

The average of %(t) = e l ' lfit / h {cc\ x )e~ liit l h is 



-10 1 2 Re (A); -1 1 



,(i)) 2 ' 2 = C £(/|n)(a x ) n , n <(n'l/K 



(E a -E m ,)t/h 



(D5) 



From Eq. (JHSJ) we have (n|/) = Xnek z g z (k z )s 2 F n (k x ) 
and the matrix element (a x ) n ^ n > is straightforward. The 
summation in (v x (t)) 2 ' 2 over si and s 2 gives two non- 
vanishing terms. We have 



dkxdk z Xnek z Nn'e'kzXn'e'k z X eters: d 

)t/h v 



-e'E, 



{v x (t)) 2 ' 2 = -c J2 

n,n' ,e,e' 

(S n ',n+1 + <V,n-l) \9z{k z )Y 

There is Xnek, = i 1 + emc 2 ) / (2E n , k J and N ne k z Xnek* = 
e/ (2E n .k,). Performing the summation over n', integra- 
tion over k x and replacing in the second term n — > n + 1 
we obtain 



(D6) 



(v*{t)f 



c r°° 

- ^2 Vn + l U n ,n+x I dk z \g z (k z )\ 2 x 

i(eE n , h!s — e'E„ + i ifej8 )t/ft 



emc 

E n _i. 



e'huj 



E n 



- Vn+l U n+hn / dk z \g z (k z )\ 2 x 

■n c d J —OO 

(££„ + i jt2 -e'£„, t Jt/t 

(D7) 



ehu 

E n ,k z J E n +i^k, 



There is 



I te l e i{eE n -eE n+1 )t/h 



(E 



n+l 



E n )t 



(E n +i + E n )t 



(D8) 



and the summations over the two terms with single e and 
e' cancel out. Rearranging terms in Eq. (|D7[) we obtain 
the same result for (v x (t)) 2,2 as in Eq. (|D2[) . Calculations 
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FIG. 9: Calculated time evolution of dynamic averages: a) 
<i!(t)> 2 ' 2 , b) (A 2 {t)) 2 '\ c) (A+(t)) 2 ' 2 , d) (A+(t)) 2 ' 2 , as given 
in Eqs. (fE2| - (|E5| . for 2+1 DE. Trap parameters as in Fig. [5] 
simulated gap frequency f2 = 2tt x 4000 Hz. Packet param- 



0.63A C , d a 
plotted for < t < 8 ms. 



0.57A C , fco* = 0.999A; 



Motion is 



for (vy(t)) ' are similar to those given above. Since a y 
has both positive and negative anti-diagonal elements, 
the expression for (v y (t)) 2 ' 2 in Eq. (|D7|) has two terms 
with opposite signs. Therefore the summation over e, e' 
cancels out the terms containing cosine functions, which 
appear in Eq. (|D8[) . and only terms with sine function 
survive. After rearranging these terms we also recover 
Eq. (|Dip . This way we showed that the average velocity 
obtained from the differentiation of (y(t)) 2 ' 2 and (x(t)) 2 ' 2 
are equal to the average values of operators (ca y (t)) 2 ' 2 
and (ca x (t)) 2 - 2 . 



Appendix E 



Below we analyze the structure of electron motion. 
Time evolution of the average values of A(t) and A + (t) is 
equivalent to the evolution of jour sub-packets: {A\{t)), 
(A 2 (t)}, (Af(t)}, (At(t)}, see Eqs. ff4fl-ffS]i. We take 
the packet = (0, f(r), 0, 0) T and follow the method 

similar to that presented in the calculation of (*4i) in 
Eq. (|T3[) . For simplicity we consider the 2+1 Dirac equa- 
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tion setting \g z (k z )\ 2 — > S(k z ), which gives 
(ii(t)) 2 ' 2 = 



E Vn + l U n . n+1 ^ e 



i(eE ni0 -E n+1 , )t/h 



1 + e' 



1 + ee 



E n +i,o 



.(El) 



Performing the summation over e, e', and writing E n 
E n ,o, = C^n+i — E n )/h, w% — {E n+ i + E n )/h, U n 
y/n + 1 U n , n+ i, and £/t = y/n + 1 U n+ i, n we obtain 



tX! W « cos («#) + T l- cos K*) - sin + iTI+ sin (w*t) } 



cos (u£t) + T+r cos (wf t) + iT++ sin «t) + iT+I sin (w^t) } 



(Mt)) 2 ' 2 


n 


{Mt)) 2 ' 2 


n 


(Af(t)f> 2 









(E2) 
(E3) 
(E4) 
(E5) 



r 



where we used the notation 

..2 



rnc" mc E„ 
si + s 2 — h s 3 — h s 4 



-En 



+ 1 



(E6) 



with si, S2, S3, S4 = ±1. Each of the terms in Eqs. (|E2[) - 
(|E5j) contains sine and cosine functions with the cyclotron 
and ZB frequencies. The structure of these terms is sig- 
nificantly different. To see this we consider the non- 
relativistic limit: E n+ \ ~ E n ~ mc 2 . Then the mo- 
tion of sub-packets (.Ai(i)) 2,2 and (^(t)) 2,2 reduces to 
the cyclotron motion, while the averages (A.2(t)) 2,2 and 
(„4^(<)) 2,2 vanish. The above sub-packets describe natu- 
ral components of the electron motion in a magnetic field. 
The direct averaging of (A(t)} or (A + (t)), as presented 
in the previous sections, allows us to calculate the evolu- 
tion of the physical quantities but it does not exhibit the 
structure of the motion. The exact operator results, as 
given in Eqs. ([55 ]) - ([38 ]) . provide a deeper understanding 
of this structure. 

In Fig.[3]we plot time evolutions of the four sub-packets 
<ii(i)> 2 ' 2 , (i 2 (t)) 2 ' 2 , (A+(t)) 2 > 2 and (A+(t)) 2 ' 2 , calcu- 
lated with the use of Eqs. (|E2[) - (|E5|) for simulated gap 
frequency f2 = 2ir x 4000 Hz. At low magnetic fields, the 
components (^(i)) 2 ' 2 and (A% (t)) 2 ' 2 are much smaller 
than (ii(i)) 2 ' 2 and (i+(t)> 2 < 2 . Note that (A^t)) 2 ' 2 
spins in the opposite direction to (Ai(t)) 2 ' 2 , and simi- 
larly for („4 2 (t)) 2,2 and (A^it)) 2 ' 2 . The four components 
of motion are persistent for the 2+1 Dirac equation. 



Appendix F 

In this Appendix we discuss the relation of our work to 
that of Barut and Thacker (BT, Ref. |7|) concerned with 



the same subject. Barut and Thacker calculated the ZB 
of relativistic electrons in the presence of a magnetic field 
at the operator level. Their work was the first treatment 
of this subject but, in our opinion, it suffered from a few 
deficiencies. 

Barut and Thacker considered the time dependence 
of electron motion introducing from the beginning its x 
and y components [in our notation, cf. Eqs. (JH) and (fT0|) 
and Appendix lA] rather than A and A + operators. This 
choice was unfortunate since A and A + satisfy separately 
important operator equations (|25[) and (|26l) . in which 
B = exp(—i£lt)A and B + = A + ex.p(+iClt) operators 
stand at the RHS and the LHS, respectively. The op- 
erators x and y do not satisfy such equations and, 'forc- 
ing' x and y to satisfy the corresponding relations, BT in- 
troduced the frequency 012 — y / '2(mc 2 ) 2 — (hu) 2 (in our 
notation). The problem here is that for hui > y/2mc 2 
this frequency becomes imaginary leading to solutions 
growing exponentially in time. In our treatment no such 
problem occurs since all the frequencies are of the form 
uj n = (E n+ i^ z ± E n< k a )/h, i.e. they are real for all mag- 
netic fields. 

The calculation of BT gave only two interband ZB 
frequencies and two intraband (cyclotron resonance) fre- 
quencies contributing to the electron motion. On the 
other hand, we obtain two series of intraband and in- 
terband frequencies because the Gaussian wave packet, 
which we use for the averaging procedure, includes nu- 
merous Landau eigenstates in a magnetic field. On the 
other hand, BT did not introduce a wave packet project- 
ing their operator results on the ground electron state. 
In contrast to our approach the procedure of Barut and 
Thacker uses the proper time formalism. 
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